/*
    Linear Programming Library
    
    Created by Team Paragon (Shena Hoffmann, Arthur Powalka, Nasir Kenarangui)
    CSCE 689
    
    Simplex.cpp
    Implements the functions declared within Simplex.h. Used for performing operations related to the simplex algorithm.
*/

#include "Simplex.h"

// ............ CONSTRUCTORS ............ //

Simplex::Simplex() { }

Simplex::Simplex(const LinProg& lp) : LinProg(lp) //dispatches lp to LinProg copy constructor
{ }

Simplex::Simplex(const Simplex& lp) : LinProg(lp)
{ }

Simplex::Simplex(bool max, vector<double> objective, vector<vector<double>> constraints, vector<string> inequalitySigns, vector<double> inequalityVals, vector<bool> constrainedVars)
{
    this->max = max;
    this->objective = objective;
    this->constraints = constraints;
    this->inequalitySigns = inequalitySigns;
    this->inequalityVals = inequalityVals;
    this->constrainedVars = constrainedVars;
    verifyLinProg();
}

// ............ DESTRUCTOR ............ //

Simplex::~Simplex()
{
}

// ............ HELPER FUNCTIONS ............ //

void Simplex::setComplexObjective(vector< complex<double> > newObjective)
{
    complexObjective = newObjective;
}	

void Simplex::findPivot(Simplex &data) 	//should return a triple (row,col,isnan)
{
	if(checkRealLessThanZero(data)==1 || checkImagLessThanZero(data)==1)	//check if this statement is necessary
	{
		int col;
		if(checkImagLessThanZero(data)==1)	//imaginary first since it represents the coefficient of M and is much larger
		{
			col=indexOfLowestImag(data);
		}
		else
		{
			col=indexOfLowestReal(data);	//case is ensured if the first doesn't occur. the top level if statement ensures this
			if (imag(data.complexObjective[col])!=0)		//if it is -a+Mb it doesn't attempt to pivot around this point. This is not a potential pivot since -a+Mb is a positive number
			{
				data.notFound=1;
			}
		}
		data.col=col;
		vector<double> tempVec(data.constraints.size());
		for(int h=0; h<data.constraints.size();h++)
		{
			tempVec[h]=real(data.inequalityVals[h])/real(data.constraints[h][col]);//A is in standard form so all entries are real... only possible none real entry is b[0] which we don't need			
		}	
		
		int row=min_element(tempVec.begin(),tempVec.end(),minPosComp())-tempVec.begin();
		if (data.constraints[row][col]<=0)	//logic error changed. see simplex.cpp for former
			data.notFound=1;
		else
			data.row=row;			
	}
	else
	{
		data.notFound=1;
	}
}

void Simplex::basisExchange(Simplex &data)		//column reduction
{
	int row=data.row;
	int col=data.col;
	//cout<<"pivot:"<<"index:["<<row<<"]["<<col<<"] value: "<<data.constraints[row][col]<< endl;
	double denom=data.constraints[row][col]; //if i don't use this then for_each doesn't work correctly... strange error where data.A[row][col] value changes ... should be the same 
	for_each(data.constraints[row].begin(),data.constraints[row].end(),[&denom](double &value){value/=denom;});
	data.inequalityVals[row]/=denom; 	//row+1 because i added an element b to represent the right side of the optimization function
				//consider instead adding an extra variable to store the right side of the objective function... that way b can be a vector of doubles instead of complex<double>'s
	for(int i=0;i<data.constraints.size();i++) //zero's out the other rows in the column
	{
		if(i!=row)
		{
			double temp=(-1)*data.constraints[i][col];	//this has to be here otherwise some columns don't change at all
			for(int j=0;j<data.constraints[i].size();j++)
			{
				data.constraints[i][j]=temp*data.constraints[row][j]+data.constraints[i][j];
				//newA(i,:)=((-1)*A(i,col)).*newA(row,:) + A(i,:);
				
			}
			//newb(i+1)=((-1)*A(i,col)).*newb(row+1) + b(i+1);
			data.inequalityVals[i]=temp*data.inequalityVals[row]+data.inequalityVals[i];
		}
	}
	double neg=-1;			//this is necessary, above something similar occurs but the compiler converts (int)*(double) into double ... (int)*(complex<double>)!= complex<double>
	complex<double> temp2=neg*data.complexObjective[col];		//this is also necessary
	//cout<<data.complexObjective[col]<<" vs "<<temp2<<endl;
	for(int j=0;j<data.complexObjective.size();j++)	//do row operation to objective function
	{
		//cout<<" temp2: "<<temp2<<"data.A[row][j]: "<<data.A[row][j]<<" data.C[j]: "<<data.C[j]<<endl;
		
		data.complexObjective[j]=temp2*data.constraints[row][j]+data.complexObjective[j];
		
	}
	data.b0=temp2*data.inequalityVals[row]+data.b0;
}	

bool Simplex::checkRealLessThanZero(Simplex &data)	//checking C
{	
	//add tolerance code to round very small numbers to 0
	bool flag=0;
	for_each(data.complexObjective.begin(), data.complexObjective.end(), [&flag](complex<double> value){if(real(value)<0) flag=1;});	
	return flag;
}

bool Simplex::checkImagLessThanZero(Simplex &data)	//checking C
{
	//add tolerance code to round very small numbers to 0
	bool flag=0;
	for_each(data.complexObjective.begin(), data.complexObjective.end(), [&flag](complex<double> value){if(imag(value)<0) flag=1;});	
	return flag;
}

int Simplex::indexOfLowestReal(Simplex &data)
{
	vector< complex<double> >::iterator index;
	index = min_element(data.complexObjective.begin(), data.complexObjective.end(),realLess());
	return index-data.complexObjective.begin();	//converts the iterator to int
}

int Simplex::indexOfLowestImag(Simplex &data)	
{
	vector< complex<double> >::iterator index;
	index = min_element(data.complexObjective.begin(), data.complexObjective.end(),imagLess());
	return index-data.complexObjective.begin();	//converts the iterator to int
}

bool Simplex::allPositiveCoefficients(Simplex &data)
{
	int flag=1; //true that all coefficients are positive
	for(int k=0;k<data.complexObjective.size();k++)	
	{
		if(imag(data.complexObjective[k])<0) 
		{	
			flag=0;	
			break;	//break loop
		} //imaginary part is just a representation not an actual imaginary number . it represents infinity 
		// this is called big M method
		else if(imag(data.complexObjective[k])==0 && real(data.complexObjective[k])<0) //if imag part is positive then the number is positive 
		{
			flag=0;
			break;
		}
	}
	return flag;	
}

//creates an identity matrix in the event objective function has initial M's
void Simplex::convertToIdentity(Simplex &data)
{
	vector<int> colsWithImag;
	int index=0;
	for_each(data.complexObjective.begin(),data.complexObjective.end(), //written to find which columns have an M in it
	[&colsWithImag, &data, &index](complex<double> value)
	{
		if (imag(value)!=0) 	
			colsWithImag.push_back(index);
		index++;
	});

	for_each(colsWithImag.begin(),colsWithImag.end(),[&data](int column)
	{
		for(int j=0;j<data.constraints.size();j++)
		{
			if(data.constraints[j][column]==1)
			{				
				double neg=-1;
				complex<double> temp2=neg*data.complexObjective[column];
				for(int i=0;i<data.constraints[j].size();i++) //can set j to zero here since all rows will have the same number of elements
				{
					data.complexObjective[i]=temp2*data.constraints[j][i]+data.complexObjective[i];
				}
				data.b0=temp2*data.inequalityVals[j]+data.b0;
			}
		}
	});
}

vector< complex<double> > Simplex::convertObjectiveToComplexForm(vector<double> oldObjective)
{
    vector< complex<double> > newComplexObjective;
    for (int i = 0 ; i <oldObjective.size(); i++)
    {
        newComplexObjective.push_back(oldObjective[i]);
    }
    return newComplexObjective;
}

// ............ CONVERSION FUNCTIONS ............ //

Simplex Simplex::convertToSimplexForm()
{
    verifyLinProg();
    Simplex StdSmx;
    StdSmx.complexObjective = StdSmx.convertObjectiveToComplexForm(objective);
    StdSmx.constraints = constraints;
    StdSmx.inequalitySigns = inequalitySigns;
    StdSmx.inequalityVals = inequalityVals;
    StdSmx.constrainedVars = constrainedVars;
    StdSmx.max=max;
    int Osize = StdSmx.complexObjective.size();
    int Csize = StdSmx.constraints.size();
    double k =0;
    complex <double> l= 0;

    // change min to max
    if (max == false)
    {
        for (int i = 0 ; i < Osize; i++)
        {
            l = StdSmx.complexObjective[i];
           StdSmx.complexObjective[i] = -l;
        }
        StdSmx.max=true;
    }
    
    // flip sign by changing inequality
    for (int i =0; i<Csize ; i++)
    {
        if ( StdSmx.inequalityVals [i]<0)
        {                
            for (int j =0; j<Osize ; j++)
            {
            k=StdSmx.constraints[i][j];
            StdSmx.constraints[i][j]= -k;
            if(StdSmx.inequalitySigns[i]=="<=")
            {
                StdSmx.inequalitySigns[i] = ">=";
            }
            else if (StdSmx.inequalitySigns[i]==">=")
            {
                StdSmx.inequalitySigns[i] = "<=";
            }
            else
                StdSmx.inequalitySigns[i] = "=";
            }                
        }
    }
    k=0;
    // add slack
        
    //new version of method to add slack variables
    for(int i=0;i<Csize;i++)	//equation i
    {
        if(StdSmx.inequalitySigns[i]=="=")
        {
            StdSmx.constraints[i].push_back (1);
            StdSmx.complexObjective.push_back({0,-1});
            for(int k=0;k<Csize;k++)
            {
                if( k!= i)
                {
                    StdSmx.constraints[k].push_back (0);		
                }
            }
        }
        else if(StdSmx.inequalitySigns[i]=="<=")
        {
            StdSmx.constraints[i].push_back(1);
            StdSmx.complexObjective.push_back(0);
            for(int k=0;k<Csize;k++)
            {
                if( k!= i)
                {
                    StdSmx.constraints[k].push_back (0);		
                }
            }
        }
        else if(StdSmx.inequalitySigns[i]==">=")
        {
            StdSmx.constraints[i].push_back(-1);
            StdSmx.constraints[i].push_back(1);
            StdSmx.complexObjective.push_back(0);
            StdSmx.complexObjective.push_back({0,-1});
            for(int k=0;k<Csize;k++)
            {
                if( k!= i)
                {
                    StdSmx.constraints[k].push_back (0);	
                    StdSmx.constraints[k].push_back (0);	
                }
            }
        }
    }

    // change to sign constrained
    for (int i = 0 ; i < Osize; i++)
    {
        if (StdSmx.constrainedVars[i]!=true)
        {	
            for (int j =0; j<Csize ; j++)
            {
                if (StdSmx.constraints[j][i]!=0)
                {
                 k = StdSmx.constraints[j][i];
                 StdSmx.constraints[j].push_back(-k);
                }
                else
                   StdSmx.constraints[j].push_back(0);
            }
            StdSmx.complexObjective.push_back(-StdSmx.complexObjective[i]);
        }
    }
    return StdSmx;
}

// ............ SOLVING FUNCTIONS ............ //

tuple<vector<double>, string, double> Simplex::solve()
{	
	Simplex data = convertToSimplexForm();
	double negative=-1;
	for_each(data.complexObjective.begin(),data.complexObjective.end(),[&negative](complex<double> &value){value*=negative;});
	data.b0=0;
	bool initImag=0;
	//printTableau(data);
	for_each(data.complexObjective.begin(),data.complexObjective.end(),[&initImag](complex<double> value)
	{if(imag(value)!=0) initImag=1;}); //checks if complexObjective has initial M's
	if (initImag==1)
	{
		convertToIdentity(data);	//converts to identity
	}
	//printTableau(data);
	findPivot(data);
	while(!data.notFound)
	{
		basisExchange(data);
		
		//printTableau(data);
		findPivot(data);
	}
	//add code here to indicate whether solution is optimal, infeasible, or unbounded
	//unbounded if the RHS of objective function has an imaginary number
	//infeasible if coefficients of objective function are not all positive
	//optimal if the coefficients of objective function are all positive and the RHS of objective function is real	

	bool checkPositive=allPositiveCoefficients(data);
	vector<double> solution;
	string feasible;
	if(checkPositive==1 && imag(data.b0)!=0)
	{
		feasible="not feasible";
	}
	else if(checkPositive==1 && imag(data.b0)==0)
	{
		feasible="feasible";
		solution=returnSolution(data);
	}
	else if(checkPositive==0)
	{
		feasible="Unbounded";
	}
	double objValue=0;

	if(feasible=="feasible")
	{	
		objValue=real(data.b0);	
		if(max==0)	//If the original problem is a minimization, our algorithm converts it to a maximum by multiplting the objective by -1
		{		//here we are multiplying it by -1 to get the solution to the original problem	
			objValue=(-1)*objValue;
		}
	}
	auto solutionData=make_tuple(solution, feasible, objValue);	
	return solutionData;
}

vector<double> Simplex::returnSolution(Simplex &data)
{
	vector<int> columns; //will indicate which variable from complexObjective indices
	vector<int> rows;    //will indicate which value from constraintVals indices	
	
	for(int r=0;r<data.complexObjective.size();r++)
	{;
		if(real(data.complexObjective[r])==0 && imag(data.complexObjective[r])==0)
		{ 
			bool isOne=0;
			for(int y=0;y<data.constraints.size();y++)
			{
				if(data.constraints[y][r]==1)
					isOne=true;
			}
			if(isOne)
				columns.push_back(r);	//stores the indices where complexObjective is 0
		}
	}	
	
	for_each(columns.begin(),columns.end(),[&data,&rows](int value)
	{
		for(int k=0;k<data.constraints.size();k++)
		{
			if(data.constraints[k][value]==1)
				rows.push_back(k);				
		}
	});
	
	vector<double> solution(data.complexObjective.size(),0);

	for(int t=0;t<columns.size();t++)
	{
		solution[columns[t]]=data.inequalityVals[rows[t]];
	}
	// in the event there are multiple solutions, the code above duplicates the solution. The code below makes the duplicate ones zero
	if (columns.size()!=data.inequalityVals.size()) //indicates multiple solutions
	{
		for(int h=0;h<columns.size();h++)
		{
			for(int j=h+1;j<columns.size();j++)
			{
				bool equalColumns=0;
				for(int g=0;g<data.constraints.size();g++)
				{
					if (data.constraints[g][h]==data.constraints[g][j])
						equalColumns=1;
					else
						equalColumns=0;
				}
				if(equalColumns)
				{
					solution[j]=0;				
				}
			}
		}
	}	
	
	return solution;
}

// ............ PRINTING FUNCTIONS ............ //

template<typename type>
void Simplex::printMatrix(vector< vector<type> > &matrix)
{
	cout<<endl;
	for(int r=0; r<matrix.size();r++)
	{
		for(int c=0; c<matrix[r].size();c++)
		{
			cout.setf(std::ios::left, std::ios::adjustfield);
			cout<<setw(23)<<matrix[r][c];
		}
		cout<<endl;
	}
	cout<<endl;		
}

template<typename type>
void Simplex::printVector(vector<type> vec)
{
	for_each(vec.begin(), vec.end(), [](type value){cout<<value<<"\t";});
	cout<<endl;
}

void Simplex::printTableau(Simplex &data) //find reference for this
{
	vector<vector<string> > matrix;
	ostringstream stream;
	vector<string> strVersion;
	for_each(data.complexObjective.begin(),data.complexObjective.end(),[&strVersion,&stream](complex<double> value){stream<<value;strVersion.push_back(stream.str());stream.str("");});
	stream<<data.b0;	
	strVersion.push_back(stream.str());
	matrix.push_back(strVersion);
	stream.str("");
	strVersion.clear();
	
	int j=0;
	for_each(data.constraints.begin(), data.constraints.end(),[&](vector<double> vec)	//adds Matrix A and its b's to the string matrix
	{
		for_each(vec.begin(),vec.end(),[&](double value){stream<<value;strVersion.push_back(stream.str());stream.str("");});
		stream<<data.inequalityVals[j++];
		strVersion.push_back(stream.str());
		matrix.push_back(strVersion);
		stream.str("");
		strVersion.clear();
	});

	printMatrix(matrix);
}

void Simplex::printProgram(Simplex p)
{
    string variable;
    if (p.max)
    {
        cout << "Maximize ";
        variable = "x_";
    }
    else
    {
        cout << "Minimize ";
        variable = "y_";
    }
    for (int i = 0; i < p.complexObjective.size(); i++)
    {
        if (i > 0) cout << " + " << p.complexObjective[i] << variable << i+1;
        else cout << p.complexObjective[i] << variable << i+1;
    }

    cout << endl << "Subject to:" << endl;

    for (int i = 0; i < p.constraints.size(); i++)
    {
        for (int j = 0; j < p.constraints[i].size(); j++)
        {
            if ( j > 0) cout << " + " << p.constraints[i][j] << variable << j+1;
            else cout << p.constraints[i][j] << variable << j+1;
        }
        cout << p.inequalitySigns[i] << p.inequalityVals[i] << endl;
    }
}

// ............ OVERLOADED OPERATORS ............ //
Simplex Simplex::operator=(const Simplex &lp)
{
    max = lp.max;
    objective = lp.objective;
    constraints = lp.constraints;
    inequalitySigns = lp.inequalitySigns;
    inequalityVals = lp.inequalityVals;
    constrainedVars = lp.constrainedVars;
}
